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Abstract 

This paper presents experimental evidence and theoretical models supporting that dry friction 
stick-slip is described by self-organized criticality. We use the data, obtained with a pin-on-disc 
tribometer set to measure lateral force to examine the variation of the friction force as a function 
of time. We study nominally flat surfaces of aluminum and steel. The probability distribution of 
force jumps follows a power law with exponents u in the range 2.2 - 5.4. The frequency power 
spectrum follows a 1//° pattern with a in the range 1 - 2.6. In addition, we present an explanation 
of these power-laws observed in the dry friction experiments based on the Robin Hood model of 
self organized criticality. We relate the values of the exponents characterizing these power laws to 
the critical exponents D an v of the Robin Hood model. Furthermore, we numerically solve the 
equation of motion of a block pulled by a spring and show that at certain spring constant values 
the motion is characterized by the same power law spectrum as in experiments. We propose a 
physical picture relating the fluctuations of the force with the microscopic geometry of the surface. 

PACS numbers: 05.65.+b,46.55.+d,64.60.Ht 



1 



I. INTRODUCTION 



There are experimental and theoretical studies suggesting that certain far from equilib- 
rium systems with many degrees of freedom naturally organize in a critical state, releasing 
energy through rapid relaxation events (avalanches) of different sizes, these sizes being dis- 
tributed according to a power law probability density. Examples of such behavior are found 
in earthquakes biological systems the stock market H, rainfall p|, and friction 

H 0, H[- All these phenomena share the features of the prototypical sand pile model || 
for which the concepts of self-organized criticality (SOC) were first proposed. Recently, the 
possibilit y o f SOC 10] in systems presenting stick-slip due to dry friction has been under 
scrutiny |lll |. In particular, Slanina |(| presented theoretical attempts to explain dry fric- 
tion in terms of SOC. However, the central question remains unanswered: to what extent 
is dry friction stick-slip a manifestation of SOC? The clarification of this issue has practi- 
cal as well as fundamental implications. From the practical point of view, the power law 
exponents could be used as parameters to characterize friction and wear of surfaces. From 
the fundamental point of view, there is a growing interest to understand systems driven far 
away from equilibrium from a single unifying principle. In addition, there is not yet a full 
understanding of the dissipation mechanisms in friction. Particularly, an overall description 
of the topography of the interface would be useful. 

In the present study, we first present experimental results on stick slip in dry friction using 
a pin-on-disc arrangement, set to measure lateral forces. The probability distributions of 
force jumps sizes and the corresponding frequency power spectra for aluminum and steel are 
examined for evidence of SOC. Second, we present a theoretical explanation of the observed 
power laws based on the Robin Hood model || [l2| , which has been successfully used in the 
past to study dislocation motion and friction. 



II. EXPERIMENT 

The pin-on-disk tribometer used for these experiments is shown in Fig. ^ This config- 
uration was chosen because it allows for easy replacement of the contacting surfaces, and 
it is the standard method for measuring friction and wear in unlubricated and lubricated 
contacts. The apparatus uses a 2.54 cm diameter disk and a spherical pin with a 0.95 cm 
radius machined on its end. The pin is attached to a load arm that is mounted on a gimbal 
supported at the center through which a load applied at the end of the arm is transferred 
to the contact zone. A strain-gauge is mounted at the end of the arm to monitor tangential 
friction force. The tangential force is monitored at a sample rate of 1,000 scans/sec and 
conversion is done using a 16-bit data acquisition card controlled by Lab VIEW. Data is 
recorded to a text file for later analysis. Frictional force measurements are done on match- 
ing aluminum and steel (M50) pin-and-disc tribometers. The signal is collected at lKHz 
during 16min, thus collecting 10 6 points. The first quarter of each data set, or about 4 
minutes, is discarded to assure that steady state is reached. In order to drive the system 
very slowly away from static equilibrium, we select slow rotational speeds in the range 10-20 
RPM. Each disk is used for up to four tests by changing the radial position of the pin on 
the disk. Loads for aluminum range from 250 g to 1000 g. Steel is studied with a 1000 g 
normal load between pin and disk. Figure 121 shows a typical tangential friction force time 
series. It shows force jumps of various sizes. We first construct the probability distribution 
of force jumps. Force jumps are taken as those events corresponding to negative changes in 
the tangential force. 

We next obtain the probability distributions of the tangential friction forces corresponding 
to aluminum under various loads and M50 steel. Results of such analysis are presented in 
Figs. H3 through |3 We observe an approximate linear behavior on the double logarithmic 
plots suggesting the power law behavior of the distributions P(F) ~ F 11 with the exponents 
fi in the range between 2.2 and 5.4. 
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We also compute the power spectra of the tangential friction force time series. We divide 
the original data (2 19 ) into 2 8 statistically independent non-overlapping data sets of 2 11 
points each. Next, we calculate 2 8 individual power spectra and finally average them to 
obtain the resulting power spectrum. Specific double-logarithmic plots are shown in Figs. [TJ- 
ITU1 The power spectra follow power laws with exponents 1.0 < a < 2.6. The results of the 
probability distribution and power spectra analysis are summarized in Table I. 



III. THEORY 

Friction is believed to occur on the atomic scale due to asperities on the surfaces in 
contacts. The simplest way to model such asperities is to use lattice models in the spirit 
of the invasion percolation [3], the sandpile model ||, Bak-Sneppen evolution model |14j . 
and Zaitsev's Robin Hood model [l2T |. In general, these models are too crude to provide 
quantitative agreement with all experimental statistical quantities. However, quantities 
such as distribution of avalanche sizes are usually described by power laws characterized 
by exponents belonging to a few distinct universality classes. These exponents may be 
compared with experimentally found ones. If several models belong to the same universality 
class, it is reasonable to study the simplest among them, since it will provide the clearest 
understanding of the physical mechanisms of the phenomenon under study. One such simple 
and elegant example is the Robin Hood model, which was originally proposed for dislocation 
movement ^2 an d later was adopted for modeling dry friction || by Slanina who added to 
the original Robin-Hood model several parameters aimed to better capture the dry friction 
mechanism, but essentially obtained the same type of behavior as the original model had. 
Here we return to the original Robin Hood model due to its simplicity. 

The model consist of a ci-dimensional lattice. Each site i on this lattice at any time step 
n is characterized by the height hi(n) which we assume to be the height of an atomic scale 
asperity at a given point of the interface between two bodies in contact. Here we present 
the model for d — 1, which is an appropriate choice to model slide friction but the analytical 
treatment is the same in any dimension, although the physically relevant cases are only 
d = 1 and d = 2. As the bodies slide against each other, the asperity with the maximal 
height is destroyed and some random number of atoms from this asperity is distributed 
among the neighboring asperities. To be specific, at each time step the site i with maximal 
height h m (n) = max hi(n) is found and the new heights are determined according to the 
following rule: h m {n + 1) = h m {n) — r{n) and h m ±i(n + 1) = h m ±i{n) + r(n)/2, where r{n) 
are independent random variables uniformly distributed between and 1. (Robin Hood 
determines the richest merchant in the market, robs him by a random amount r(n) and 
distributes it equally among the neighbors without leaving anything for himself). If we 
assume periodic boundary conditions so that the sites with i = and i = L are equivalent, 
the total amount of matter J2i=o^i( n ) * s conserved and we can assume it to be zero. The 
distance between the surfaces at a given site % can be determined as h m {n) — hi{n). 

The particular details of the model such as the distribution of r(n) or the rule of dividing 
it among neighbors can vary, but the model still retains its SOC behavior . The critical 
exponents appear to be sensitive to the details of dividing r(n), for example, an exactly 
solvable asymmetric model (in which all the profit is given to the site on the left) [l5[ 
belongs to a different universality class. 

It has been shown [l6l. [lTj that in a wide class of depinning SOC models, all the critical 
exponents can be expressed in terms of the two main exponents : avalanche dimension D and 
correlation exponent v. The avalanche of threshold ho is defined as sequence of time steps 
during which the hight of the maximal asperity is above h . Namely, if h m (no) > h and 
h m (no + s) > h while for n < n < n + s, h m {n) < h , the sequence n = n + 1, ..,n + s 
is called a punctuating avalanche of threshold ho and mass s. The avalanche dimension 
describes how the avalanche mass s, scales with the horizontal dimension of the avalanche 
R. To be more precise, the mass distribution of forward avalanches with threshold ho, scales 
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as 

P.{s) ~ s- T >g 8 (s(ho - h c ) Dv ) (1) 
and the distribution of the avalanche horizontal size, R, scales as 

P R (R)~R- T «g R (R(h -h c y), (2) 

where h c « 0.114 is the critical height, 

r s = l + (d-l/u)/D (3) 

and 

TR = l + d - 1/u (4) 

are Fisher exponents first introduced to characterized cluster distributions in percolation 
theory [lal. w hile g s and g R are exponentially decreasing cutoff functions. It has been 
suggested ]19j . that the Robin Hood model belongs to the same universality class as the 
linear interface model, for which the values (D = 2.23 and r s — 1.13 in d — 1; D — 2.725 
and t s = 1.29 in d = 2) are given in Ref . |17l| . Using these values and Eq. (jSj), one gets 
i/ = 1.41 for d = 1 and z/ = 0.83 for d = 2. 

It can be shown that after the initial equilibration number of time steps T ~ L D , any 
initial shape of the interface hi(0) reaches a steady state such that very few N(L) "rich" 
sites have hi(n) > h c ~ 0.114, where 

N{L) ~ (5) 

and 

d f = d-l/u (6) 

plays the role of fractal dimension of rich sites. Only those few rich sites have a chance to 
be robbed. The chance P m (h m ) that at a given time step, the maximal height is equal to 
h m decreases for an infinite system [l6[ as 

P m (h m ) = (h m - h c y~\ (7) 

where the exponent 

7 = -tZ) (8) 

characterizes the dependence of the average avalanche size on its threshold h : (s) ~ (ho — 

The distribution of heights of the poor sites converges to a smooth distribution on the 
interval [h c — 1, h c ], while the distribution of the rich sites converges to the distribution with 
a power law singularity 

P h (h)~(h-h c )- dv . (9) 

This result is not presented in Refs. [HI [lTj but can be justified by the following heuristic 
arguments. Indeed, the number of sites with h > h scales as the number of the active sites 
in an avalanche of threshold h , and thus scales as R df (ho), where R(h ) is the cutoff of the 
avalanche distribution (J2J) which scales as 

R(ho) ~ (h - h c y v . (10) 

Thus the probability that h > ho scales as {ho — h c )~ d f u and the probability density of h = ho 
scales as 

Ph(ho) ~ (^o - he)-*"- 1 = (h - h c y dv . (11) 
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We can assume that the friction force, F(n), at a given time step is proportional to the 
number Ph(h m (nj)Ah of asperities with heights between h m (n) — Ah and h m (n): 

F(n) = F x P h {h m {n))Ah, (12) 

where Ah is the interaction distance of atomic forces acting between the two surfaces and 
F\ is a proportionality coefficient, corresponding to the surfaces interaction force at the 
asperity. Accordingly, the distribution of the friction forces P(F) satisfies the equation 
P(F)dF = P m (h m )dh m , where the random variables F and h m are linked by Eq. (fI2j). 
Taking into account Eqs. Q and (JE) we have dh m /dF ~ dF~ 1/du /dF ~ i?-V^-i. Finally 
Eq. © yields 

rlh 

P(F) = P m (h m (F))^ ~ F -to-W*> F -V*>-i = (13) 

where 

fi = (D + l/v)/d. (14) 

For d = 1, using values of Ref. we have -P(-F) = -F M with /j = 2.94 which is consistent with 
the experimental observations of the density of jump sizes presented here and in Ref. [20]. 
For d = 2 we have \i = 1.96. 

In order to test this theoretical predictions, we perform simulations of the one dimensional 
Robin Hood model. Starting at n = with a flat interface hi(0) = 0, and selecting the first 
site to rob at random, after T steps we get all L sites of the interface updated at least 
once. Measuring the average (T) for many independent runs for different system sizes, and 
plotting it versus L in a double-logarithmic scale (Fig. HTjl . we can obtain the avalanche 
dimension D as the limit of the successive slopes of this graph for L — > oo. 

A typical shape of the interface at time n > T is presented in Fig. One can see that 
the height of the majority of sites do not exceed the critical value h c ~ 0.114. Interestingly, 
the majority of rich sites with heights above the critical barrier are localized in the vicinity 
of the richest site. 

Figure EH shows the histogram of all the interface heights Ph{h) collected over many time 
steps after the system has reached the steady state and the histogram of the heights of the 
robbed sites P m (h m ). One can see that while Ph(h) dramatically increases as h — > h£, no 
sites below the critical value are robbed and P m {h m ) — > as h m — > h+. In order to find 
the exponents governing the behavior of these distributions near the critical point, we plot 
these quantities in a double logarithmic scale as functions of h — h c (Fig. fl3*b). 

Finally we determine the time series of forces, defined as the number of heights between 
h m (n) and h m (n) — Ah as function of time. (Fig. 1X1)1 . The histogram of this time series is 
presented in Fig. [J5]in a double logarithmic scale. The slope of this plot is // = 3.0, which 
is consistent with the theoretical prediction ()14|) . 

Note that the time series F(n) is slightly correlated, which can be demonstrated by 
the negative slope of its power spectrum Sp(f) ~ f~ a in the log- log scale (Fig US)). The 
explanation of this phenomenon is based on the fact that values h m (n) fluctuate in the 
vicinity of h c in a non-trivial way, so that h m (n) become less than h c + e at time steps 
n separated by intervals distributed according to Eq. (JTJ). This is because these intervals 
coincide with avalanches for the threshold h Q = h c + e. The values of h m (n) below h c + e 
correspond to the large values forces F(n) and thus the intervals between the forces F(n) 
above certain threshold are also distributed according to Eq. (JTJ). It can be shown |17L l2ll|tIiat 
the exponent a of a time series generated by peaks separated by intervals of zero signal 
distributed according to a power law as in Eq. ((TJ) is equal to r s — 1 for 1 < r s < 2. Thus 
according to Eq. (J5J) a = (d — l/u)/D ~ 0.13. Indeed, the numerical data of Fig. ITfil give 
a ~ 0.14 in a very good agreement to the above theoretical prediction. However, this value 
of spectral exponent is much smaller than the values observed experimentally which are in 
the range between 1 and 2.6. 
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This difference is to be expected since the materials in contact as well as the strain gauge 
have finite elastic constants and inertia which produce a time delay between the applied 
force and the displacement record by the tribometer and lead to an effective integration 
of the input force time series. We would expect that if the materials were infinitely stiff 
then the experimental force power spectrum should agree with the theoretical predictions. 
Therefore, we construct a mechanical model of a tribometer, that accounts for these effects. 

We assume that the pin of the tribometer contacts the sample at time t at a point with 
coordinate x(t) and it is dragged along the sample by the strain gauge spring with spring 
constant k attached to the body of the instrument moving along the sample with constant 
velocity Vq, which is equivalent to the rotational speed of the disk. The force measured by 
the tribometer is thus k[vot — x(t)], which fluctuates as the pin moves against the sample 
with velocity v(t) = dx/dt and acceleration a(t) = d 2 x/dt 2 . The equation of motion of the 
pin is thus 

ma = (vot — x)k — F(t, v), (15) 

where F(t,v) is the friction force generated by the highest asperity of the sample and m is 
the mass of the pin. 

Now our goal is to relate F(t,v) with the input from the Robin Hood model. Note, that 
the physical time t is not directly proportional to the time step n of the Robin Hood model, 
but is equal to the sum of the durations of each time step 

n 

* = !>, ( 16 ) 

i=l 

where the durations £, are the times needed for the pin to travel a characteristic distance 
Ax, which is the linear size of each asperity. We assume that if the pin moves along the 
sample by Ax , the current asperity is destroyed and the landscape of the contact between 
the pin and the sample is rearranged according to the rules of the model. Thus the time step 
n t of the Robin Hood model, corresponding to a given moment of time t can be determined 
as nt = mt[x(t)/ Ax], where int[...] denotes the integer part of the expression in the brackets. 

If v = 0, F(t,v) = sign(t> t — x) mm[bF(n t ), k(v t — x)], where F{n t ) is the input from 
the Robin Hood model, and b is some material and load-dependent constant. If v ^ 0, 
F(x,v) = sign(t> )[bF(n t ) + rjv], where r\ > is some dissipative constant. Constant b is 
proportional to the load and depends on the elastic properties of the material. Introducing 
dimensionless variables by x' = xj Ax and t' = tv^/Ax, we arrive to a dimensionless equation 

a! = (t' -x')k' -F'(t',v'), (17) 

where k! = kAx 2 jmv 2 . and F'(t',v') is the same as F(t,v) but the constants b and r] are 
changed by V = bAx/mv 2 , rf = rjAx/mvo. 

Thus, there are three independent dimensionless parameters of the model: k! , b' and 
i]'. Varying these parameters, we found a wide region in the parameter space in which the 
power spectrum of the model resembles the experimental one. A typical example of the 
spectrum for k! = 0.001, b' = 0.3 and r/ = 0.01 is shown on Fig. (jl7j) . The frequency of the 

resonance peak is determined by y/k' /2tt ~ 5 • 10~ 3 . The peak becomes more pronounced as 
we decrease rf. The increase in rj' also increases the slope of the spectrum. The increase of 
b' at given k' increases the frequency region in which the power spectrum follows the power 
law but it also increases the absolute value of the slope closer to 2, a characteristic value 
of the Brownian motion. In general, an integration of the time series corresponds to the 
increase of the spectral exponent by 2, so the integration of the white noise produces the 
Brownian noise. The observed spectral exponent a = 1.45 suggests that in a certain range 
of parameters, our model acts as the fractional integrator of the input signal. For a very 
stiff spring and large dissipation (k = 1, b = 0.1, r] = 1) the output signal of our equation is 
not much different from the input time series F(n t ) and we recover the small value of the 
spectral exponent a = 0.14. 
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IV. CONCLUSIONS 



We present experimental results and theoretical arguments that support the presence of 
self organized criticality in dry sliding friction. The experiments are pin-on-disk friction 
force traces of aluminum-aluminum and steel-steel systems. In both cases and for a variety 
of normal loads, the distribution of the friction force jumps and the frequency power spectra 
are power laws. The theoretical arguments are based on the application of the Robin Hood 
model to the friction problem. This model provides rules by which the surface profile changes 
as a function of time. The model introduces a height h that we interpret physically as the 
height of the asperity. At each time step, atoms from the highest asperity are distributed 
among neighboring sites. We use the known distribution of heights and of maximal heights 
h m of the Robin Hood model to obtain the time series of the friction forces created by the 
asperities. As the maximum height fluctuates near the critical value, the number of smaller 
asperities whose heights are within the reach of atomic forces also fluctuates, diverging as 
h m comes close to the critical value h c . These smaller asperities correspond to the contact 
sites and are responsible for the friction force. Specifically, the friction force is proportional 
to the number of contacts. Thus we propose that the friction force at a given time step is 
proportional to the probability density of the interface heights at the current value of the 
maximal height. The statistical distribution of the friction forces is studied both numerically 
and analytically. 

We also find that the large forces are bunched in time. This is due to fluctuations of 
the maximal heights above a constant critical height. When the maximum heights return 
to the critical value, the forces become large. Thus, the surface waxes and wanes between 
a situation of large force due to many asperities acting, and a situation of smaller force in 
which only few asperities are in contact. 

In addition, we use the time series of forces as an input to the Newton's equation which 
describes the kinematics of the pin. For stiff or massless materials, the experimental distri- 
bution of force jumps should coincide with the theoretical distribution of forces. However, 
materials have finite mass and elasticity and thus the experimentally measured friction forces 
differ from the actual forces at the contact. To investigate these effects, we solve this equa- 
tion numerically for different values of the parameters and find good agreement with the 
experiment. 
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FIG. 1: Pin-on-disk tribometer. The arrow points at the location where the spherical pin and the 
disk touch. The disk lies horizontally while the pin attached to the arm, rests above it. 



TABLE I: The values of exponents n characterizing the power law behavior of the distribution of 
the force jump sizes and spectral exponents a characterizing the power spectrum of the friction 



force time series for different materials and loads. 

Material Load fx a 

__ i, ),),)„.. ^ 

Al 250g. 5.4 1.0 

Al 750g. 2.2 1.5 

Al lOOOg. 3.2 1.3 
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FIG. 2: Typical signal from the tribometer. The effective spring constant of the apparatus is 
lg//xm, giving the largest force jumps as a few hundred /j,m. 
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FIG. 3: Probability density for jump size distribution on steel M50. 
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FIG. 4: Probability density for jump size distribution on aluminum with a normal load of 250g 
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FIG. 5: As in Fig. |I] with normal load of 750g. 
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FIG. 6: As in Fig. 0] with normal load of lOOOg. 
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FIG. 8: Power spectrum of jump sizes for an aluminum sample with a normal load of 250g. 
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FIG. 9: As in Fig. El with a normal load of 750g. 
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FIG. 11: Double logarithmic plot of the average equilibration time (T) versus system size L = 
2 3 ,2 4 , ...,2 12 . The inset shows the successive slopes of the main graph versus 1/L. The intercept 
D = 2.23, agrees with the data of Ref.pjj- 
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FIG. 12: A typical shape of the steady state interface. The horizontal line shows the critical 
height. 
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FIG. 13: (a) The semi-logarithmic plot of the distribution of heights Ph{h) of all sites (solid 
line) and the distribution of heights -P m (/i m ) of robbed sites (dashed bold line). Vertical dotted 
line shows the position of the critical height h c = 0.114. (b) Double logarithmic plot of the same 
quantities plotted as functions of h — h c . The slopes of the curves in the fitted regions are in 
agreement with Eq.@ -dv = -1.4 and Eq. (J7J) v(D - 1) = 1.72. 
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FIG. 16: Power spectrum of the time series F(t) presented in Fig. (|14j) . The slope of the straight 
line fit is —a — 0.14. 
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